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An efficient real space method is derived for the evaluation of the Madelung's potential of ionic 
crystals. The proposed method is an extension of the Evjen's method. It takes advantage of a general 
analysis for the potential convergence in real space. Indeed, we show that the series convergence is 
exponential as a function of the number of annulled multipolar momenta in the unit cell. The method 
proposed in this work reaches such an exponential convergence rate. Its efficiency is comparable to 
the Ewald's method, however unlike the latter, it uses only simple algebraic functions. 



PACS numbers: 

I. INTRODUCTION 

Since 90 years a number of methods have been pro- 
posed to calculate the electrostatic potential in ionic crys- 
tals. These methods can be separated into two categories, 
the direct summation methods and the indirect summa- 
tion ones. The former uses a real space summation of 
the electrostatic potential generated by the ions within a 
finite volume (£). However, when enlarging the volume 
£ , such partial summations are conditionally convergent. 
The convergence depends on the specific shape of £. fn 
addition, when achieved, the convergence is quite slow. 
The indirect summation methods do not present these 
drawbacks since the long range part of the potential is 
calculated in the reciprocal space. Indeed, the summa- 
tion is divided into two parts, a short range one, eval- 
uated by a direct summation in real space and a long 
range one evaluated in the reciprocal space. Among 
these methods, the most widely used is the Ewald's 
method!, which is actually considered as the reference 
for Madelung potential calculations. 

Despite its quality the Ewald's method is not easily 
usable in different domains of physics. This is for in- 
stance the case in clusters ab-initio calculations used for 
the treatment of strongly correlated systems, the study of 
diluted defects in materials or adsorbates or for QM/MM 
type of calculations. For this types of calculations, real 
space direct summation methods are used. There is thus 
a need for efficient and accurate techniques for the deter- 
mination of the Madelung potential in real space. 

The convergence problems found in real space summa- 
tion are linked to the shape of the summation volume, 
£ , and more specifically the charges at its surface. In 
order to insure the convergence of the summation, the 
surface charges are renormalizcd. Several methods have 
been proposed for this purpose. 

The most common and simple one is the Evjen's 
method^. This method uses a volume £, built from a 
finite number of crystal unit cells, and renormalizes the 
surface charges by a factor 1/2, 1/4 or 1/8 according 
whether the charge belong to a face, edge or corner of £. 
This method insures, in most cases, the convergence of 



the electrostatic potential when £ increases. However, in 
some cases such as the famous CsCl it does not converge 
to the proper valued 

Other authorsM proposed to renormalize not only the 
surface charges, but also the charges included in a thin 
skin volume. The adjustment of the rcnormalization fac- 
tors are, in this case, numerically determined so that to 
reproduce the exact potential, previously computed us- 
ing the Ewald's method at a chosen set of positions. Such 
a method presents the advantage of reaching a very good 
precision. However several drawbacks can be pointed 
out : i) the previous calculation of the electrostatic po- 
tential using the Ewald's method at a large number of 
space positions, ii) the necessity to invert a large linear 
system to determine the rcnormalization factors and iii) 
finally the fact that the latter are not chosen on physi- 
cal criteria. Indeed, this last point induces the possibility 
that the renormalization factors can be either larger than 
one or negative. It results that even if the electrostatic 
potential is very accurate at the chosen reference posi- 
tions, its spatial variations can be unphysical and thus, 
the precision can strongly vary when leaving the reference 
points. 

Marathe et a£ suggested a physical criterion, based on 
the analysis of the convergence of the real space summa- 
tion, for the choice of the renormalization factors. In- 
deed, it is known that the direct space summation con- 
verges to the proper limit if the volume £ presents null 
dipole and quadrupole momenta^. The authors of ref- 
erence [H showed, on the simple example of a linear al- 
ternated chain, that it is possible to find a finite number 
of charge renormalization factors allowing the cancella- 
tion of these two multipolar momenta. They also assert 
that the cancellation of additional multipolar momenta 
increases the speed of convergence. Unfortunately they 
did not prove this affirmation and more importantly, they 
did not proposed a practical way to determine the renor- 
malization factors in order to reach this goal. 

In the present paper we propose a systematic method 
for the determination of the renormalization factors al- 
lowing the cancellation of a given number of multipo- 
lar momenta as well as a careful analysis of the direct 
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space summation convergence as a function of the num- 
ber of canceled multipolar momenta. The next section 
will present the convergence proof, section 3 will develop 
the method for determination of the renormalization fac- 
tors and section 4 will present the optimization of the 
method and illustration on a typical example. 



II. CONVERGENCE ANALYSIS 
A. Potential at a point 

As already mentioned in the introduction, several pa- 
pers already exist on this subject. However the results 
are only partial and there is not complete analysis of the 
convergence issue. We will thus present in this section a 
global analysis of the electrostatic potential convergence 
in a real space approach and an estimation of the error. 

We want to evaluate the limit of the following series 



V(f ) = lim V n (r ) = lim V V 
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where R is a vector of the Bravais's lattice, j refers to 
a charge qj located at the position Fj of the unit cell C. 
{£ n , n £ N} is a set of volumes such that 
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For the sake of simplicity we require that the set of £„ 
also presents the following conditions 
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The well known problem of this series is that the limit 
depends on the particular choice of the unit cell C and of 
the volumes £ n . However, it has been demonstrated that 
this conditional convergence disappears if one considers 
a unit cell with zero dipolar and quadrupolar momenta-. 
We will consider in the following that the cell C fulfills 
this condition. 

The error AV n (fo) on the electrostatic potential eval- 
uation, V n (fo), can be written as 
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For large values of n, AU„(fb) can be evaluated by a 
multipolar expansion. For practical reasons, we will use 
an expansion expressed in spherical coordinates. Indeed, 
for a given order, this expansion contains less terms than 
the usual multipolar expansion based on Cartesian co- 
ordinates. The use of the later multipolar expansion is 
still possible but is more complex (see refQ). In spherical 



coordinates, the multipolar expansion of the error made 
on V n (r ) reads^ : 



AK (ro) = E E E M lm (ro)^^ (6) 

R££„ 1 m=-i 



where M.i m (fo) are the multipolar momenta of unit cell 
C at r . They can be expressed as 

M lm (f ) = J2 <liVj - ro\ l Yr m (dj,^) (7) 
jec 

(R, 9, 4>) arc the spherical coordinates of the Bravais's 
vector R and (pj , 8j , <f>j ) the spherical coordinates of 
fj — r*o- The Yf m are Schmidt semi- normalized spherical 
harmonics : 



where P™ are Legendre functions. 

In order to overvalue the error, one needs to overvalue 
the spherical harmonics. We thus consider the addition 
formula : 



Y l m (6,<l>)Y l - m (9',<f>') = P l (co S i) (9) 

771= — I 

where P; is a Legendre Polynomial and 7 is the angle 
between (R, 9, <j>) and (R 1 , 9', 0'). It comes for 9 = 9' and 



\Yr(e,^\ 2 =Pi(i) = i 
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and thus 
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(ii) 



Using this result, one obtains the following overvalue 
for the momenta : 



\Mi m \ < Q(r + a) 1 



(12) 



where a is the typical size of C (i.e. the diameter of 
the circumsphere of C) and Q is the sum of the absolute 
values of its charges 



£1*1 



(13) 



One should notice at this point that Y l m (9, </>) has the 
same parity as I. The contributions of R and —R unit 
cells to AV n (ro) thus cancel when I is odd. Using equa- 
tions [TT] and [T21 one gets the following overvaluation : 



\AV n (f )\<Qj2(4k + l)(r + a) 
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where 2p is the order of the first even, non-zero momen- 
tum. 

Let us now overvalue the sum over the 1/R powers 
by a volume integral. For each cell C(R) located at R, 
the norm of the position vectors r belonging to the vol- 
ume C(R) is smaller than R + a. One can thus overvalue 
1/R 2k+1 by 
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lo being the volume of the unit cell C. If TZ n is the radius 
of the insphere of £ n , it comes 
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and 
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k—p 



The later sum converges for TZ n large enough, i.e. for 
TZn > »"o + 2a. The decreasing function (4k + l)/(2k — 
2) can be overvalued by its value in k = p. Further 



summation over k leads to the following expression 
|AF„ (f ) \<a p fnD 5 ( ^ )2 \ )2p -x («) 
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The electrostatic potential at ro thus converges as l/TZ l r ~ 2 
where / 
cell C. 



where / is the first, even, non-zero momentum of the unit 



B. Difference of potential between two points 

In several applications, as for instance in cluster ab ini- 
tio calculation, the problem depends on the spatial vari- 
ations of the potential and not on its absolute value. In 
such cases it is sufficient to cancel the dipolar momentum 
of the unit cell in order to ensure the convergence of the 
calculation. The convergence rate can also be expected 
to be faster than for the calculation of the potential at a 
point as we will show in this section. 

Let us overvalue the error made on the calculation of 
a difference of potential between two points located at 
f Q + fi and f - fi : 
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where (R-, 0-, </>_) and (R + ,8 + ,4> + ) are the spherical 
coordinates of R — r\ and R + f i respectively. As in the 
preceding section (R, 8, <fi) will be the spherical coordi- 
nates of R and (ri , 0i , 0i ) those of fi . 

In order to express the previous expression as a func- 
tion of 1/-R, we use following expansion of solid spheri- 
cal harmonics (for simple derivation see ref. [l(| see also 
ref.EIIll) : 

P^(cos 0_) e !Mfa r n fL + n- m\ 

R L+1 - 2^ R L+n+l L — M J 

x (-l) M - m P^ +n (cos0) e im * 
x P r f -™(cosfli) e^ M ~ m ^ (21) 

where the sum over m spans all integer values. Neverthe- 
less, only a finite number of terms will contribute, since 
P\ n = if \m\ > I. Setting I = L + n and introducing 



spherical harmonics leads to : 
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Considering R + in this equation instead of i?_ is equiv- 
alent to the transformation 

0i — > 7T-01 
4>1 > 7T + 01 

that results in an overall (— 1)'~ L factor. Inserting rela- 
tion [221 into eq. [20] and inverting the summation over I 
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and L leads to the expansion : 
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Considering the parity of Yi m , one can see from eq. [53] 
that the contributions from cells located at R and — R 
cancel when I is odd. Moreover, due to the last term of 
cq. [Ml the Ai m coefficients are zero when I — L is even. 
Only terms with I even and L odd have a non zero con- 
tribution, thus only momenta Mlm^) with odd order 
will contribute to the error. The consequence is that the 
first non-zero contribution in equation 1231 corresponds to 
I = 2p + 2 where 2p + 1 is the first, non-zero, odd mo- 
mentum of the unit cell. 

Let us now find an overvalue of the Ai m terms. It is 
easy to show using a recurrence relation on the values of 
m and M, that if \m\ < I, \M\ < L and \M -m\ <l-L, 
the following relation holds : 



I — m \ / / + m \ (I 
L-Ml \L + m) ~ \L 



(25) 



Using the previous overvaluation of the momenta (eq. [T 
one obtains : 
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As in previous section, the sum over 1Z can be overval- 
ued by a volume integral (cf. eq[H]). Overvaluation of 
the error thus reads : 
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The sum over k converges if lZ n is larger than ro + 2a + 
r\. It can be calculated using derivative of power series. 
After simplification, one obtains : 



\AV n {r ,ri)\ < b p 



TZ 2 n (r + a + r 1 ) 2 P+ 2 

2p 



(TZ n -r -2a- r x ) 2 {H n - a) 



where 



8ttQ 16p 2 + 40p + 25 
bp = x 

uj 2p 



(28) 
(29) 



As one increase the size of the set of charges, the dif- 
ference of electrostatic potential between two points con- 
verges like l/7?4 -1 , where I is now the first, odd, non-zero 
momentum of the unit cell C. This convergence is slightly 
faster than the convergence of the absolute value of po- 
tential. When the order of the first non zero momentum 
is even, the convergence rates differ by a factor 1/i? 2 , 
otherwise they arc similar. 



III. PARTIAL CHARGES 

As depicted in the previous section, convergence can 
be considerably increased if one cancels several multi- 
polar momenta of the unit cell. In general the Evjcn 
method allows to only cancel the dipolar momentum, and 
thus provides a convergence of the potential differences 
in In order to really take advantage of the former 

property, one needs a method allowing the cancellation 
of several multipolar momenta. 

In this section we will establish a method to con- 
struct unit cells with a chosen number of zero multipolar 
momenta. The method, based on the usage of partial 
charges, is general and can be applied to any Bravais's 
crystal. 

Let (a ,b ,c ) be the lattice vectors of the Bravais's 
crystal, and Cq the associated unit cell. In order to in- 
troduce partial charges, we consider a larger cell C; of 
dimensions (I x ao,l x bo, I x Co), that we will refer as 
the 11 construction ceW . The construction cell thus con- 
tains I 3 original unit cells Cq which positions in Ci can be 
labeled by p, q and r indices, ranging from 1 to I. 

If we note n c the number of charges qi in the orig- 
inal cell Co, the cell Ci now contains n c x I 3 charges. 
These charges will be corrected by a factor X pqr (where 
i refers to the charge qi). When on rebuild the lattice 
using the construction cells Ci, the cells overlap, and the 
final charge at position fj corresponds to the superposi- 
tion of partial charges from several construction cells. It 
is straightforward to show that the condition to retrieve 
the nominal value of the charges qi reads : 



E 



(1 < i < n c 



(30) 



At this stage, considering the latter n c equations, the 
cell Ci contains n c {l 3 — 1) free parameters that could be 
used to cancel multipolar momenta. For the sake of sim- 
plicity and generality (i.e. for the method not to depend 
on the particularity of a given crystal), we will impose 
further conditions on the X pqr coefficients. 

We first reduce the problem to a one dimensional prob- 
lem by setting : 



Xp' 1 X q 



b.i 



(31) 



where the three coefficients A°' 1 , X q ' 1 and Xp 1 are used to 
cancel multipolar momenta of the ID problems obtained 
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when the cell Co is respectively projected on the three 
axes of the crystal. For each one dimensional problem, 
the construction cell contains I projected unit cells. The 
condition on the coefficients now reads : 



P =i 



(u> = a, b, c) 



(32) 



It is easy to show that, if these coefficients cancel a fixed 
number of multipolar momenta in each one dimensional 
problems, the X pqr coefficients will cancel the momenta 
of same order in the original three dimensional problem. 

We further impose to the coefficients to only depend 
on the fractional coordinates (a,-, A'j 7i) of the charges g,-, 
in the corresponding direction: 



a:, 



Ap(aj) \ q (/3i) A r (7i) 



(33) 



The Xp(x) functions are thus the same for the three di- 
rections and for all charges. As a consequence their ex- 
pression is the same for all crystals. For a given value of 
x, and considering the condition for the reconstruction of 
the crystal (eq. [32]) , we are left with I — 1 degrees of free- 
dom. We thus impose to the X p (x) functions to cancel 
I — 1 multipolar momenta. This can be done by setting 
the momenta created at the center of the construction 
cell by a unique charge q : 



qui Y^ X p( x ) [% + P- 1 - 2 J = 1 u o m Lk (34) 

(0 < k < I - 1) 



where uq = ao,&o> or cq, x is the fractional coordinate 
of the charge and m;,fc are constant values (independant 
of the crystal specifications). The equation obtained for 
k = corresponds to the condition for the reconstruction 
of the crystal (mi t o = 1). The momenta of the construc- 
tion cell, can thus be obtained by summing the contri- 
butions of all charges. As the unit cell is neutral, these 
contributions cancel out. 

The equations [33J and [3J] thus define sets of partial 
charges that allow to construct cells with I — 1 zero mul- 
tipolar momenta. The shape of these partial charges de- 
pends on the choice of the mi,k constants values. In order 
to find the most reasonable choice of partial charges, we 
search for m^k constants that satisfy the following phys- 
ical conditions : 

1. Vp, X p (x) e [0,1]. 

2. The partial charges vary continuously, i.e. 
Vp, X p (x) are continuous and A p (l) = A p +i(0). 

3. Vp, X p (x) decreases monotonously when moving 
away from the center of the construction cell. 

4. The values of the mi k are as small as possible. 



The latter condition ensures that the partial charges are 
larger in the center of the construction cell and smaller on 
its edges, and hence that this cell is close to the original 
one. 

We did not find a way to derive the solution of this 
problem in a general way, for any value of I. We thus 
determined the solution for fixed values of / up to I = 6. 
In all these cases the first X p function presents the same 
shape : 



J-i 



Xi(x) 



(35) 



We reasonably assume that this expression is valid for 
any value of I. As we will see, it is possible to show, 
a posteriori, that the A p functions fulfill the first three 
conditions. 

We will now determine the function X p for any I, using 
the above expression of Ai and relation [3J] In equa- 
tion [331 the rtik,i constants can be replaced by the value 
of the momenta obtained for x = : 



^Ap(x) (x+p — 1 
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E- 



A p (0) 
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(36) 

with < k < I — 1. The matrix of this linear system is 
the transpose of a Vandermonde Matrix. Inversion of the 
system leads to : 



A p (x) 
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q=l 



a 9 (o) n 



x + i — q 
-P 



(37) 



The A p (0) coefficients can now be obtained from the 
expression of Ai. Let us consider the latter equation eval- 
uated for p = 1 , and for I integer values of x : 



q— 1 v 



(1 < j < l) (38) 



Replacing Ai(n) by its expression and inverting this lin- 
ear system leads to the expression of A g (0) coefficients : 



9-1 



A 9 (0) = E(-!) J 



3=1 



l\(q-J-l) 



(1-1)! 



(39) 



Finally, using eq. 1371 and 1391 one obtains the expression 
of the X p (x) which are polynomial functions of order / — 1 
in x (see fig. [J) : 

3=1 j=l V/ V ' i=l. i^p y 

(40) 

These functions are segments of the uniform sum distri- 
bution , i.e. the distribution Pi(x) of the sum of / uniform 
variates on the interval [0, 1] : 

X p {x) = Pi{x+p- 1) (0 < x < 1 and 1 < p < I - 1) 

(41) 
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FIG. 1: (color online) X P {x) functions obtained for values of I ranging between 2 and 10. Intervals correspond to cells Co that 
compose the construction cell C;. In each intervals, fractional coordinates x range between and 1. 




x+p-1/2-1/2 

FIG. 2: (color online) Renormalization of the charges at the "left" edge of a crystal fragment, obtained for values of I ranging 
between 2 and 10. The origin of the abscises corresponds to the position of the edge obtained when using the original cell Co 
(dotted line). 



From this relation, it is obvious that the X p (x) functions 
satisfy the first three conditions mentioned above. 

We now consider a fragment of crystal made of n con- 
struction cells Ci . As these cells are composed of Z original 
cells Co, they partially overlap, and the size of the frag- 
ment corresponds to n + I — 1 cells Co- n 1 = n — I + 1 
cells Co in the center contains charges with the nominal 
values qi , and I — 1 cells Co on each side of the fragment 
contains partial charges. The latter partial charges are 
proportional to the coefficients : 

p 

fi p (x)=J2\(x) (l<P<i~l) (42) 

9=1 

These coefficients are represented on fig. [3] The abscise 
values have been shifted by (Z — l)/2, so that the origin 
corresponds to the position of the effective edge of the 
fragment (i.e. the position of the edge obtained when us- 
ing n original cells Co without partial charges). One can 
see that the renormalization of the charges is relatively 
small. Indeed, even in the case I = 10, this renormal- 
ization is weaker then 1% for charges at distances larger 
than 2uq, (uo = ao> bo or Co). 



IV. OPTIMIZED METHOD 

We first use the cells C; to illustrate the convergence 
of the standard, real-space calculation of the potential 
established in section [TTJ In order to observe the general 
behavior of the different methods, we chose the a quartz 
structure, that possesses a reasonable number of atoms 
per unit cell and not too many symmetries. A cell C; is 
constructed for a fixed number of p. The cells are used to 
produce sets of charges of increasing size. The potential 
is calculated at the position of the twelve atoms of the 
central unit cell Co. 

Fig.[3]a) represents the maximum of the error made on 
these potential values as a function of the width of the set 
of charges, and fig. [3] b) represents the maximum of the 
error made on the sixty six differences of potential. As 
expected the error decreases as a power function of the 
width D of the set of charges. It fully agrees with eq. [T51 
and 1251 In particular, the fact that the cancellation of a 
momentum of odd order do not increase the convergence 
rate of the potential at a point, clearly appears. Similarly, 
the cancellation of even order momentum do not improve 
the convergence rate of the calculation of differences of 
potential. 

One can see from the previous figures that this stan- 
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FIG. 3: (color online) a) Maximum of error made on the 
calculation of the potential at the position of the atoms of 
the central unit cell. D — n' + 21 — 2 is the width of the 
sets of charges in each crystallographic direction, where n' 
is the width of the central zone where the charges have their 
nominal values. I is the order of the first non zero momenta of 
the cell Ci used to construct the set of charges, b) Maximum 
of error made on the calculation of the difference of potential. 



FIG. 4: (color online) a) Maximum of error made on the 
calculation of the potential at the position of the atoms of 
the central unit cell. The abscise corresponds to the width 
D = n! +21 — 2 of the sets of charge in each crystallographic di- 
rection, where n' is the width of the central zone with nominal 
charge values, b) Maximum of error made on the calculation 
of the difference of potential. 



dard approach is not the more efficient. Indeed the in- 
crease of the number of zero multipolar momenta clearly 
yield a faster convergence rate than the increase of the 
volume of the system for a fixed value of I. Let us there- 
fore fix the width n' of the volume containing the nominal 
charge, and let I increase. The variation of the maxi- 
mum error made on the potential and on the the poten- 
tial differences are respectively represented on fig. |4] a) 
and|l]b). One sees that the present approach leads to an 
exponential convergence of the potential in both cases. 
The convergence is very fast, since an increase of the set 
of charges width by two unit cells results in a precision 
increase by a factor better than ~ 10. Increasing the 
number n' of central cells without partial charges has a 
small influence on the convergence speed. For n' > 3 
the method even becomes less efficient since increasing 
n' increases the size of the total set of charges. The best 
convergence is obtained for n' = 3. It corresponds to the 
case where the cell in which the potential is calculated is 
surrounded by one shell of cells with the nominal charge 
values. 



Finally we compare our method to the famous Ewald's 
method which mixes calculation in real space and recip- 
rocal space. This method introduces Gaussian distribu- 
tion of charge exp(— a 2 r 2 ), where the a coefficient can 
be adjusted. Increasing a coefficient increases the con- 
vergence rate of the real space sum, but slows down the 
sum in reciprocal space. A width of Gaussian propor- 
tional to the characteristic length of the unit cell, which 
corresponds to an = cj -1 / 3 , is generally assumed to give 
a good compromise. We calculated the error made on the 
value of the potential using the Ewald's method for dif- 
ferent values of a around an ■ The results are represented 
on figure as well as the error of our method obtained 
for n' = 3. Figure [5] reports the error on the potential 
as a function of the number of construction cells used in 
the calculation. Let us point out that, while this variable 
is pertinent for the global convergence rate analysis, for 
each charge, the Ewald's method requires an error func- 
tion evaluation resulting in a non negligible pre- factor, 
not present in our method and not taken into account 
in figure [5] One sees that the convergence rate of the 
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FIG. 5: (color online) Error on the Madelung's potential 
evaluation using the present method (bold black solid curve) 
and the Ewald's method with different parameters a (colored 
dashed curves). 

present method is comparable with the Ewald's method. 
If one is only interested in the potential evaluation at 
a single point, the Ewald's method with an optimal a 
parameter is somewhat faster than the present one. One 
the other hand, once the renormalization have been com- 
puted, the value of the potential at any other point of the 
of the central area can be calculated with a similar pre- 
cision at little cost. More important, properties using 



potential integrals or complex potential functions can be 
more easily evaluated since our method used only alge- 
braic functions. 



V. CONCLUSION 

Number of authors have searched for a fast converging 
method for the evaluation of the electrostatic potential 
in real space. Similarly, many works where done yield- 
ing partial results on the convergence rate of such real 
series. The present work fills the gaps and proposes a 
general analysis of both the convergence of the poten- 
tial at one point and of the convergence of differences of 
potential. Indeed, we gave a general and rigorous proof 
of the relation (claimed by other authors) between the 
power law convergence of the series and the number of 
zero multipolar momenta of the crystal construction cell. 

Based on these convergence analyses we derived a gen- 
eral real space method with an exponential convergence 
rate, comparable with the Ewald's method. The expo- 
nential convergence is reached as a function of the num- 
ber of canceled multipolar momenta in the construction 
cell. The crystal is indeed constructed using overlapping 
construction cells with rcnormalizcd charges. We derived 
a general analytical expression of the renormalization fac- 
tors, for any given number of zero multipolar momenta. 
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